%%
clc;clear
syms p(n) z
assume(n>=0 & in(n,'integer'))
% f = p(n+2) - p(n+1) - p(n)
f = n+1
%%
fZT = ztrans(f,n,z)

%%
syms pZT
fZT = subs(fZT,ztrans(p(n),n,z),pZT)

%%
pZT = solve(fZT,pZT)

%%
% pSol = iztrans(pZT,z,n);
pSol = iztrans(fZT,z,n);
pSol = simplify(pSol)

%%
pSol = subs(pSol,[p(0) p(1)],[1 2])

%%
nValues = 1:10;
pSolValues = subs(pSol,n,nValues);
pSolValues = double(pSolValues);
pSolValues = real(pSolValues);
stem(nValues,pSolValues)
title('Rabbit Population')
xlabel('Years (n)')
ylabel('Population p(n)')
grid on